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Abstract 

We present molecular dynamics simulations of the homogeneous 
(mechanical) melting transition of a bcc metal, vanadium. We study 
both the nominally perfect crystal as well as one that includes point 
defects. According to the Born criterion, a solid cannot be expanded 
above a critical volume, at which a 'rigidity catastrophe' occurs. This 
catastrophe is caused by the vanishing of the elastic shear modulus. 
We found that this critical volume is independent of the route by which 
it is reached whether by heating the crystal, or by adding interstitials 
at a constant temperature which expand the lattice. Overall, these 
results are similar to what was found previously for an fee metal, 
copper. The simulations establish a phase diagram of the mechanical 
melting temperature as a function of the concentration of interstitials. 
Our results show that the Born model of melting applies to bcc metals 
in both the nominally perfect state and in the case where point defects 
are present. 



1 Introduction 

Over the years, several theories explaining the mechanism of melting have 
been proposed. |2l Ej This research has by now evolved to a state where 
a clear distinction exists between two possible scenarios for the melting 
transition: a first scenario of homogeneous, or mechanical melting resulting 
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from lattice instability [Sllllini and/or a spontaneous generation of thermal 
defects, [71 IHl El Uni El and a second scenario of heterogeneous, or thermody- 
namic melting which begins at extrinsic defects such as a free surface or an 
internal interface (grain boundaries, voids, etc). [12111211111 113 CHI Through- 
out this paper we will use the term mechanical melting to describe the former 
case, which we consider here. In particular, we take the view proposed by 
Born that at the melting point a 'rigidity catastrophe' is caused by the van- 
ishing of one of the elastic shear moduh,[31 Hj C44, or C = (Cn — Ci2)/2. 
In other words, the crystal melts once it loses its ability to resist shear. 
This condition determines the mechanical melting temperature, T^, of a per- 
fectly homogeneous bulk crystal and was confirmed in extensive studies of 
fee metals. [131121120111111101 

TallonHj pointed out that a mechanical instability arises when the solid 
expands up to a critical specific volume which is close to that of the liquid 
phase (melt). In the study by Wang et al. fHlCSl of the mechanical melting 
transition of an fee solid under external stress, it was found that volume 
expansion is the underlying cause of lattice instability. Kanigel et al. [TT] 
confirmed this scenario in a simulation of fee copper in the presence of point 
defects. They showed that the critical volume at which a crystal of copper 
melts is independent of the path through phase space by which it is reached, 
whether by heating of the perfect crystal or by adding point defects to expand 
the solid at a constant temperature [TT]. 

Solids can undergo mechanical melting only if they have no extended 
defects [0], a situation which is conveniently realized in three-dimensional 
computer simulations by means of periodic boundary conditions which elim- 
inate the surface. Simulations of atomic dynamics for solids with an fee 
[THl HOI 1201 HH UOj or diamond [T7j structures show, among other things, the 
onset of a shear instability of the solid at a temperature Tg, which can exceed 
the thermodynamic melting temperature Tm by some 20%, depending on the 
details of the potential. 

Given the considerable degree of understanding of the melting process of 
fee crystals, it is of substantial interest to see if the scenario of mechanical 
melting also applies to solids having a different lattice structure. We there- 
fore decided to study mechanical melting of a bee metal, vanadium, by means 
of computer simulations. We present details of the calculations in Sec. II. In 
Sec. Ill we describe the results of the simulation of some physical properties 
of a vanadium with and without point defects. In Sec. IV we present molec- 
ular dynamics simulation results for mechanical melting in the presence of 
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point defects. Finally, in Sec. IV, we discuss the implication of our results 
for the question of the microscopic mechanism of melting. 

2 Simulation details 

We model the bulk melting transition of vanadium using the molecular dy- 
namics simulation [21] technique. The choice of vanadium has no special 
significance as we are only interested in the generic features of metallic solids 
with a bcc lattice symmetry. While various many-body potentials for fee met- 
als [nn] have been developed and thoroughly tested in numerous simulations, 
the situation with such potentials for bcc metals is not as good. This can be 
explained by the more complicated nature of the bcc metals in comparison 
with fee ones which manifests itself in the wide range of elastic constants 
(and even in the negative values of the Cauchy condition C12 — C44 < for 
some bcc metals). The packing density of atoms in a bcc lattice is smaller 
than in a fee lattice (there are 8 nearest-neighbors in a bcc lattice and 12 
nearest-neighbors in a fee). However, the second nearest-neighbor distance 
in the bcc structure is larger than the first nearest-neighbor distance by only 
about 15%. Therefore, the interaction between the second order and the first 
order nearest-neighbors in bcc metals is not negligible, even with screening. 

In addition, band structure effects are crucial for bcc metals. A simple 
approximation which assumes that the electron density can be considered as 
a superposition of atomic orbitals is successful for fee metals, but less appro- 
priate for bcc metals. Therefore, for metals with the bcc structure electron 
density is chosen to be an adjustable function, rather than a superposition 
of atomic orbitals. Furthermore, angle dependent interactions could be very 
important in bcc solids. 

For our simulations, we chose the many-body interaction potential devel- 
oped by Finnis and Sinclair [21] (FS), and modified by Rebonato et al. [22] FS 
proposed a way to incorporate the delocalized physical nature of the metallic 
bonding and the essential band character of bcc metals in a simple model. 
The FS potential involves two short ranged potentials, a cohesive one and 
a repulsive one. The cohesive potential is summed over neighbors and the 
square root of the result describes the bonding energy. The cohesive energy 
is therefore proportional to ^/z, where z is the atomic coordination number. 
The square root form is used assuming that the band energy is the sum of 
occupied one-electron levels, and according to the tight-binding model the co- 
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hesive energy is proportional to the hopping integrals between the d-orbitals. 
The repulsive potential is summed in the usual way to describe the repulsive 
core-core interactions. 

Our molecular dynamics (MD) simulations with the FS potential were 
performed using the Parinello-Rahman method and the Nose-Hoover 
thermostat. j2Sl I2n] This ensemble is identified as an isothermal-isotensional 
ensemble (NtT), j^Tj which allows simulation of fiuctuations in the shape and 
volume of the sample (here, is the number of atoms, T the temperature, 
V the volume and t is the sample tension). The shape and volume thus 
obtained were used for calculation of the shear modulus in a canonical (NVT) 
ensemble. The shear elastic moduli were calculated using the fiuctuations of 
the stress tensor. |2H] 

The samples used for the simulations contained 2000 atoms, initially ar- 
ranged as a perfect bcc crystal of size 10x10x20 unit cells. Periodic bound- 
ary conditions were applied in all three directions. Point defects were in- 
troduced either by the insertion of extra atoms between the lattice sites 
(self-interstitials) or by the removal of atoms from the lattice (vacancies). 
Newton's equations of motion were solved using a predictor-corrector algo- 
rithm. 1211 ESi Throughout this study, interactive visualization with the AViz 
program [30] was implemented. 

3 Validation of the potential and order pa- 
rameter 

To learn about the capatibility of the potential, we examined some physical 
properties of a perfect crystal. First, we calculated the thermal expansion 
at zero external pressure. We found the thermal expansion coefficient at 
low temperatures to be = (18 ± 6) x 10~^ in reasonable agreement 

with the experimental value measured at room temperature ctexp = 8.6 x 
10^^ K^^. Next, the thermodynamic melting temperature for our potential 
was calculated, using the method of Lutshko et al. [33], to be = 2220 ± 
10 K. This value is close to the experimental value = 2183 K, despite 
the fact that the FS potential was constructed by fitting its parameters to 
room temperature values of various physical properties of vanadium (lattice 
constant, cohesion energy, shear elastic moduli, vacancy formation energy, 
etc). 
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In order to test the algorithm we calculated the shear moduli as a function 
of temperature. The shear elastic coefficients decrease with temperature as 
shown in Fig.^ The accuracy of the simulations was estimated by monitoring 
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Figure 1: Variation of C (triangles) and C44 (squares) with temperature. 
The error bars represent statistical uncertainty. 

the convergence of the shear elastic moduli calculated along symmetrically 
equivalent directions. We found the difference to be approximately 10%. 

Following the validation that our potential can indeed reproduce the phys- 
ical properties of a perfect crystal with acceptable accuracy, point defects 
were introduced. These point defects are distributed homogeneously through- 
out the bulk of the solid. Only one type of point defects, e.g. vacancies or 
self-interstitials were used in each run to avoid their mutual annihilation. 

The configurations of atoms in the vicinity of a point defect inside the 
bulk at low temperatures was investigated by means of the simulated tem- 
pering method. The most energetically favored configuration of an 
interstitial was found to be the < Oil > dumb-bell split - interstitial (See 
Fig. El) with a formation energy of Ef = 4.18 ± 0.02 eV. This formation 
energy is in agreement with that of previous simulations 

To investigate the temperature dependence of the crystalline order, we 
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Figure 2: The most energetically favored configuration is found to be the 
< Oil > dumb-bell split - interstitial 

define the structure order parameter rj: 




where k = {0, 0, 27r/a} is a vector of the reciprocal lattice, ft is the position of 
atom i, N is the number of the atoms in the sample, and the angular brackets 
stands for ensemble average. For an ideal-crystal lattice at zero temperature, 
rj equals unity, while in the hquid state, r] fluctuates near zero. 

We calculated 6ri/6C, the change of the order parameter upon the intro- 
duction of point defects. Here, C is the concentration of point defects, given 
in % of the number of atoms. Fig. El shows the result of this calculation for 
small C ,and at different temperatures. The introduction of self - interstitials 
results in noticeable decrease of the structure order parameter (from rj ~ 0.6 
to ?7 ~ 0.4), while the influence of vacancies is relatively weaker. With in- 
creasing temperature, the order parameter becomes increasingly sensitive to 
the introduction of point defects, as evidenced by the increase of the absolute 
value of \Srj/6C\ with temperature. We believe that this increased sensitivity 
results from the increase of the amplitude of thermal vibration of the atoms 
in the immediate vicinity of the point defect. 
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Figure 3: Influence of point defects: vacancies (diamonds) and self- 
interstitials (squares) on the structure order parameter r] as function of tem- 
perature. Tlie concentration of point defects, C, is given as a percent of the 
total number of atoms. The error bars represent statistical uncertainty. 



Introduction of self - interstitials increases the volume of the sample (See 
Fig. while vacancies decrease the volume, as shown in Fig. The specific 
volume of point defects at various temperatures was estimated using the 
linear dependence of the volume on the number of defects, apparent in Figs. 0] 
and El The volume of the sample at a small number of vacancies can thus 
be written as 

V={N- N,a)v + N,aV,a (2) 

here, is the number of atoms in the sample, N^a is the number of vacancies, 
V is the volume per atom in a perfect crystal of vanadium and v^a is the 
volume per vacancy. A similar relation for self-interstitials can be written as 



V = Nv + NsiVs 



(3) 



where, Nsi is the number of self-interstitials and Vsi is the volume per inter- 
stitial. It is interesting to point out that the linear dependence appears to 
hold even at temperatures close to Tg. This may indicate that the concept of 
a point defect remains meaningful even under these conditions. The specific 
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Figure 4: Total volume of the sample as a function of concentration of 
interstitials at several temperatures: T = 2300 K (diamonds), T = 2200 K 
(squares) and T = 2000 K (crosses). The concentration of point defects is 
given as a percent of the total number of atoms. The error bars represent 
statistical uncertainty. 

volume of a point defect (in atomic volume units) is shown as a function of 
temperature in Fig. IHl It is seen that at temperatures above 2000K these 
specific volumes change rapidly. To a large degree, this increase can be ac- 
counted for by the rapid decrease of the elastic coefficients of the crystal in 
this temperature range. 

4 The bulk melting transition 

The prime goal of our simulations is the investigation of the role of point 
defects in mechanical melting. In the simulations of mechanical melting of 
fee metals jHl ^1 E] it was established that the key parameter controlling 
melting is the volume of the crystal. When the Born criterion is applied 
to a superheated crystal lattice it establishes the existence of a critical vol- 
ume above which the crystal becomes mechanically unstable and therefore 
undergoes a phase transformation to the liquid state or some other crystal 
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Figure 5: Total volume of the sample as a function of concentration of 
vacancies at several temperatures: T = 2300 K (diamonds), T = 2200 K 
(squares) and T = 2000 K (crosses). The concentration of point defects is 
given in (%)of the total number of atoms. The error bars represent statistical 
uncertainty. 

structure. The critical volume is coupled with a maximum superheating tem- 
perature, Tg. Simulations with fee metals P ITHl ITT] showed that this critical 
volume, Vs, can be attained by expansion caused either by heating the crystal, 
or by doping it with point defects at a constant temperature which expand 
the crystal, [TTj or by pure mechanical dilatation at zero temperature. [HI HE! 
In this sense the mechanical melting process appears to be universal, i.e. it 
is determined only by the sample expansion up to the critical volume. 

In order to verify whether the same scenario holds in the case of a bcc 
metal we carried out simulations using samples with various concentrations 
of self-interstitials or alternatively, vacancies. The initial temperature of 
each sample was chosen far below the melting point of a perfect sample, 
T ~ 0.7 Tg. As the samples were heated by gradually increasing the temper- 
ature, at some point we observed an abrupt decrease of the structure order 
parameter (see Fig. [7j), together with a simultaneous increase of the total 
energy and volume (see Fig. IH)). This event determines the mechanical 
melting temperature. The melting temperature of a sample without point 
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Figure 6: The ratio of specific volume of point defects to the specific vol- 
ume of an atom as a function of temperature: self-interstitials (squares) and 
vacancies (triangles). The error bars represent statistical uncertainty. 
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Figure 7: Typical time dependence of the order parameter during mechanical 
melting. This particular sample contained 0.25% interstitials at temperature 
T = 2475K. 
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Figure 8: Typical jump of the sample volume during the mechanical melting 
transition. This particular sample contained 0.25% interstitials at tempera- 
ture T = 2A7bK. 

defects is found to be Tg = 2500 ±20 K. Since MD simulations are plagued by 
statistical fluctuations in the temperature and volume, in practice it is very 
difficult to reach the maximum superheating temperature, Tg. Therefore, the 
accuracy in the determination of in this way is about ~ 1%. 

The same temperature, Tg = 2500 ± 12 K, was also found from a least- 
squares fit to the temperature dependence of C as shown in Fig. Q It is 
the temperature where C goes to zero. This indicates that as is the case 
for fee metals, homogeneous melting of the bcc metal results from a shear 
elastic instability. This particular value of T, applies to a crystal of vanadium 
containing no defects, and is about 280 K higher than the thermodynamic 
melting point = 2220 ± 15 K obtained for our model using the method 
proposed by J. F. Luthsko et al. [Mj 

Once point defects are introduced, it is found that becomes a function 
of their concentration. Results of simulations performed at different temper- 
atures and defect concentrations are summarized in our phase diagram (See 
Fig. El). The fact that point defects lower the melting temperature has been 
confirmed experimentally ( 7-irradiation lowers the melting point of pure 
metals by an amount proportional to the dose, and thus to the number of 
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generated point defects. [?, E]) The lowering of can be explained as fol- 
lows: Introduction of self - interstitials leads to a significant local distortion 
of the bcc lattice, and expands the volume of the solid, see Fig. ^| There- 
fore, a solid containing self - interstitials reaches its critical volume already 
at a lower temperature ( the melting temperature is lower). In contrast, 
the effect of vacancies is rather minor, at least if their concentration is small 
enough. The same effect of lowering of the bulk melting temperature induced 
by interstitials was obtained by A. Kanigel et al. JTP for copper (fee lattice). 
However, at higher concentrations of point defects the decrease of T, can 
not be explained simply by volume expansion. This is especially notable in 
the case of vacancies which decrease the volume, but at high enough con- 
centrations also lower the melting temperature (See Fig. Ej). We refer here 
to the region in Fig. El where the concentration of point defects approaches 
1%. These values are unrealistically large in comparison with the typical 
laboratory values ~ 0.001%. At these high concentrations, the concept of 
a single point defect is unclear and one should perhaps consider clusters, or 
extended defects. According to Jin et al. [20] extended defects can act as 
nucleation centers for melting. Taking this point of view, the lowering of 
with defect concentration may result from the combined effect of (a) volume 
expansion, and (b) introduction of nucleation centers for melting. Finally, it 
should be noted that the calculated phase diagram is qualitative, because of 
the finite sample size and limited simulation time. 

Our results are broadly consistent with models of defect-induced melting 
proposed by Fecht [?] and Granato. [7^ According to Fecht [?] melting is 
driven by the incorporation of point defects into the lattice. Point defects 
increase the probability of heterophase fluctuations of liquid-like clusters in 
the defective crystal and lower the Gibbs energy of the crystalline state. 
Therefore, the melting temperature decreases as the concentration of point 
defects increases. 

The configuration of point defects (self-interstitials) in a fee metals was 
exploited by Granato 7J to construct a model giving the thermodynamic 
properties of the crystalline and hquid states in a unified way. He found a 
large diaelastic softening of the shear modulus with increasing defect con- 
centration. This leads to a softening of the formation energy of interstitials, 
which, together with the large entropy contribution from the interstitialcy 
resonance modes, lowers the melting temperature. In the above discussion, 
we have emphasized the role of lattice instability in establishing a maximum 
superheating temperature at zero external pressure. However, due to ther- 
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Figure 9: The influence of interstitials (squares) and vacancies (diamonds) 
on tlie melting temperature of vanadium under periodic boundary conditions. 



mal expansion, any temperature change is accompanied by a simultaneous 
change of the volume. To decouple these two effects, we plot the dependence 
of the shear modulus C on the specific volume in Fig. ^| As this figure 
shows, the dependence of C on the specific volume appears to be universal, 
in the sense that the value of C is the same whether the volume at which 
it is calculated was reached at by thermal expansion or by insertion of point 
defects. In other words the main effect of interstitials is to expand the lattice. 
Using the data plotted in Fig.Elone can extract the value of the critical vol- 
ume, Vs{Ts)., at which the system melts homogeneously under the conditions 
of zero external stress. 

Using this method we find Vg = 14.87 ± 0.06 A^/aXom and the melting 
temperature for various concentrations of point defects. The critical vol- 
ume is close to the specific volume of liquid vanadium at the thermodynamic 
melting temperature Vug = 15.3 ±0.05 and to the experimental value [23] 
of vug = 15.2 A^. 

Similar results were obtained for copper in MD simulations by J. Wang et 
al. fSJ and by A. Kanigel et al. ^] It was found that the shear modulus 
vanishes at a lattice strain of a/a^ = 1.024, where a lattice parameter at 
Tjn = 13507^, and is the lattice parameter of copper at Tq = 300K. The 
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Figure 10: Plot of the shear modulus C against specific volume at vari- 
ous concentrations of interstitials: squares - crystal without impurities (only 
thermal expansion), diamonds - 0.05% concentration of interstitials, circles - 
0.1%, triangles - 0.15%, crosses - 0.2%. 

specific vohmie ratio of copper is (a/ao)'^ = 1.07 which is quite close to the 
value obtained for vanadium v{T^)/v{Tq) = 1.06 ± 0.01. It is natural to 
ask whether the ratio a/ao is universal, independent of lattice structure. To 
answer this question in a definitive manner, it would be useful to make similar 
simulations on other bcc metals. 

5 Summary and conclusions 

In our simulations we observed that each shear elastic modulus is a continuous 
and apparently universal function of the specific volume. The solid lattice 
can be expanded either by thermal expansion or by the presence of self- 
interstitials. The value of C at any particular volume is independent of 
the way by which this volume was reached, either by thermal expansion 
alone, or by any combination of thermal expansion and of expansion due 
to interstitials introduced into the sample at a constant temperature. The 
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elastic energy of the lattice increases until a critical specific volume Vg (close 
to the specific volume of the melt) is reached where the shear modulus C 
vanishes, triggering mechanical melting. Upon melting, the solid transforms 
isothermally and discontinuously (see Figs. [7| and |H1). 

The process that triggers mechanical melting could be similar to the one 
observed by Jin et al. [201 in simulations of the melting of a surface-free 
Lennard- Jones crystal. There, melting occurs when the superheated crys- 
tal spontaneously generates a sufficiently large number of extended defects 
(clusters of spatially correlated destabilized particles which satisfy the Born 
criterion). Those extended defects play the role of surfaces in initiating the 
melting. In our simulations, point defects, especially in large concentration 
where clusters of defects should be formed, could act as nucleation centers 
for these extended defects (molten regions) inside the solid. 

This paper was devoted to a simulation of the melting process of a homo- 
geneous bcc metal, and its comparison with a similar process in fee metals. 
It is of great interest to extend these simulations to heterogeneous melting 
which involves nucleation of the liquid phase at some preferred sites of the 
solid, for example at the free surface. This study is the subject of a forth- 
coming paper. 
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